The fraction of condensed counterions around a charged rod: 
Comparison of Poisson-Boltzmann theory and computer simulations 






o 

oo 



o 



Markus Deserno^, Christian Holm^ and Sylvio May^ 

^ Max- Planck- Institut fiir Polymerforschung, Ackermannweg 10, 55128 Mainz, Germany 

^ Institut fiir Biochemie und Biophysik, Friedrich-Schiller-Universitdt Jena, Philosophenweg 12, 07743 Jena, Germany 

(February 1, 2008) 

We investigate the phenomenon of counterion condensation in a solution of highly charged rigid 
polyelectrolytes within the cell model. A method is proposed which - based on the charge distribution 
function - identifies both the fraction of condensed ions and the radial extension of the condensed 
layer. Within salt-free Poisson-Boltzmann (PB) theory it reproduces the well known fraction 1 — 1/^ 
of condensed ions for a Manning parameter ^ > 1. Furthermore, it predicts a weak salt dependence 
of this fraction and a breakdown of the concept of counterion condensation in the high salt limit. 
We complement our theoretical investigations with molecular dynamics simulations of a cell-like 
model, which constantly yield a stronger condensation than predicted by PB theory. While the 
agreement between theory and simulation is excellent in the monovalent, weakly charged case, it 
deteriorates with increasing electrostatic interaction strength and, in particular, increasing valence. 
For instance, at a high concentration of divalent salt and large ^ our computer simulations predict 
charge oscillations, which mean-field theory is unable to reproduce. 
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I. INTRODUCTION 

Strongly charged linear polyelectrolytes use-, their 
counterions to reduce their line charge densityJd This 
phenomenon Jias led to the concept of counterion 
condensation,aB and although it was introduced a long 
time ago, varyijig viewpoints about this subject persist in 
the literature Jj'u Here we will investigalfi-its appearance 
within the commonly used cell modelWB: an infinitely 
long charged rod enclosed in a cylindrical cell together 
with its counterions — with and without added salt. 

In the salt-free case this model can be solved an- 
alytica]J|}t| within nonlinear Poisson-Boltzmann (PB) 
theory,LfEl which thus affords a particularly clear view on 
Manning-Oosawa counterion condensatian, as has been 
demonstrated by Zimm and Le BretJj'H Going beyond 
salt-free PB theory, many questions arise: How closely 
condensed and tightly bound is the "condensed layer"? 
What distinguishes condensed from uncondensed coun- 
terions? When is the mean-field level PB theory a good 
approximation for real systems? And how does the pres- 
ence of salt affect the condensation phenomenon? 

Recently Manning proposed the idea that there exists 
a clear distinction between (i) a condensed layeii^nd (ii) 
a distant, more diffuse "Debye-Hiickel" cloud.til In the 
integrated radial counterion distribution function this is 
supposed to be detectable as an inflection point, which 
separates the two regions. Here we argue that this is not 
quite in accord with PB theory without added salt: There 
is an inflection point in the distribution function, but it 
is not related to the condensation phenomenon. If how- 
ever the distribution function is plotted against logarith- 
mical radial distance, an inflection point appears which 
exactly divides the counterions into condensed and up-i 
condensed ones, as previously pointed out by Belloni.llj 
Since this feature has largely gone unnoticed in the study 
of polyelectrolytes, we found it worthwhile to present its 



derivation within PB theory in a short and explicit form 
and also to point out its practical usability. To this end, 
we present computer simulations, compare them to PB 
theory, and demonstrate that this criterion can indeed 
be extended to quantify counterion condensation even 
beyond the scope of PB theory. 

In the case of added salt the analytical treatment, 
the cylindrical PB equation is much more involved,ll3u 
so any potential inflection point criterion is more difficult 
to analyze. Still, the PB equation can be solved numeri- 
cally, suggesting that the presence of monovalent salt de- 
creases the extension of the condensed layer, but leaves 
the amount of condensed counterions largely unaffected. 
We show that this behavior is well reproduced in com- 
puter simulations for monovalent salt. We also present 
a simple criterion for determining the salt concentration 
above which the concept of Manning condensation is no 
longer meaningful. In the regime of strong electrostatics, 
high valence and much added salt the PB predictions de- 
viate qualitatively from simulational results. In particu- 
lar, the simulation shows a pronounced overcharging and 
charge oscillations, which are absent on the mean-field 
PB level. 

This paper is structured as follows: In sections || and 



[II we recapitulate the main ingredients of the PB solu- 



tion for the salt-free cell model and illustrate the connec- 
tion of Manning condensation with the infiection point 
criterion mentioned above. This is followed - in sec- 
tion IV ~ by a comparison of the salt-free PB results with 
computer simulations. In sections M and VI we discuss 
the concept of counterion condensation in the presence 
of salt and to this end derive a PB equation for the en- 
semble of constant number of salt molecules. Its results 



are compared with simulations in section VII. Details of 



our simulation method can be found in the appendix. 



II. PB-THEORY FOR A CHARGED ROD 
WITHOUT ADDED SALT 

Consider an infinitely long cylinder of radius ro and 
line charge density A > 0, which is coaxially enclosed in 
a cylindrical cell of radius R. Global charge neutrality of 
the system is ensured by adding an appropriate amount 
of oppositely charged (monovalent) counterions. 

Within PB theory these counterions are replaced by 
a cylindrically symmetric counterion density n{r) (r is 
the radial coordinate) which gives rise to an electrostatic 
potential <i>(r) satisfying the Poisson equation 
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(1) 



with e being the dielectric constant outside the cylinderll3 
and e the (positive) unit of charge. Conversely, this po- 
tential is supposed to influence the counterion density via 
the Boltzmann factor: 



n{r) — n{R) exp {/3e$(r)} 



(2) 



with the inverse temperature (3 = l/k^T and fee being 
Boltzmann's constant. Thus, the chosen zero of the po- 
tential is <P{R) = 0. 

In the following it is advantageous to change variables: 
The Bjerrum length £b — /3e^/47re provides a convenient 
scale for quantifying electrostatic interactions (it is the 
distance, at which the Coulomb energy of two elemen- 
tary charges equals /cbT), while the dimensionless Man- 
ning parameter ^ = XIb/s measures the line charge den- 
sity of the rod (it is equal to the number of elementary 
charges per Bjerrum length) . We shall only be interested 
in the strongly charged case ^ > 1. Insertion of eq 
into eq |l| results in the nonlinear PB equation which in 
terms of the reduced (dimensionless) electrostatic poten- 
tial y{r) = /3e$(r) and a screening length 1/k > with 



AT:iBn{R) reads 
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The appropriate boundary conditions for solving the PB 
equation arise from applying Gaufi' law at rg and R: 



y'(ro) - -2C/ro 



y'{R) = 



(4) 



Eq y has an analytical solution which can be written in 
the following way: 



2/(r) = -2 1n^-^cos(7ln-^)|. (5) 



nr 



The boundary conditions (^ yield two coupled, transcen- 
dental equations for the two integration constants 7 and 
Rm'- 



7 m -- — = arctan 

.Km 7 

, R 1 

7 m — — = arctan — 
.Km 7 



(6) 

(7) 



Subtracting (||) from ^ eliminates Rm and provides an 
equation from which 7 can be obtained numerically. The 



second integration constant Rm, which we will refer to 
as the Manning radius, is then given by either of these 
equations. Note also that k and 7 are connected via 
k'^R^ = 2 (1-1-7^), thus ensuring the chosen normalization 
of the potential. 

The Manning radius Rm depends monotonically on ^ 
and for ^ > 1 one finds Rm > ro. As discussed in the 
next section this is the regime in which counterion con- 
densation occurs. If ^ = 1 then Rm — ro, i.e., the Man- 
ning radius is located at the surface of the rod. A fur- 
ther decrease in ^ shifts Rm inside the cylinder and for 
^ = ln(i?/ro)/(l + ln(i?/ro)) both the Manning radius 
and 7 vanish. Even smaller values of ^ render the inte- 
gration constant 7 complex. Still, the solution (0) can be 
extended by analytic continuation over C 

Using eqs 0, ra and ra the total charge per unit length, 
Q(r), found within a cylinder of radius r £ [ro] R] can be 
determined by integration: 



Q{r)/X = l — — / dr2TTren{r) 
A Jro 

1 - (1 - ^) - ^ tan 



(7ln^ 
V -Km 



(8) 



Since n{r) > 0, Q{r) decreases monotonically from 
Q{ro) = A to Q{R) = 0. The latter follows from eq 
and is a consequence of global charge neutrality. It is 
instructive to use the quantity 



P{r) = 1 - Q{r)/X, 



(9) 



which is the probability of finding a mobile ion within 
the distance r. In other words, it is the fraction of coun- 
terions found within a cylinder of radius r. Inparticular, 
at r = Rm the last term in Q, as given in eq@, vanishes, 
giving a fraction 1 — 1/^ of ions within Rm- It can easily 
be verified that generalizing eqs m and M for counteri- 
ons with valence v reduces to replacing ^ — > ^w. Within 
PB-theory changing valence or electrostatic interaction 
strength affects the charge distribution function in the 
same way. 



III. COUNTERION CONDENSATION: 
DEFINITION AND IDENTIFICATION 

For ^ > 1 eqs y and imply the inequalities 

TT TT 

> 7 > 
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(10) 



Since the two bounds become equal in the limit R —f 00, 
they provide an asymptotic solution for 7 and give rise 
to various limiting laws, which illuminate the behavior of 
the solution in the dilute limit. In particular, the reduced 
potential for ^ > 1 becomesjla 



yW-y(^o)--2 1n--2 1n|l + (e-l)ln-|, 



which is (up to a logarithmic correction) identical to the 
potential of a rod with charge density of e/^B, i-e.. Man- 
ning parameter ^ = 1. This can be attributed to a con- 
densation of counterions onto the rod, which renormalize 



the line charge density. Indeed, the contact density n{ro) 
converges to the nonzero value 



lini ri(ro) 



A (g-1) 



(12) 



suggesting the existence of a close layer which cannot be 
diluted away. 

In order to establish an effective Manning parameter of 
1, a fraction /^ = 1 — 1/^ of all counterions would have to 
condense onto the rod. In fact, /^ is a critical threshold 
in the following sense: For 0<q;<1,^>1 and the 
radius r^ defined as r^ — r^ exp{a/(^ — 1)(1 — a)} one 
can show by using eqs g-jol that 

(13) 



lim P{ra) = af^. 



Hence, in the limit of infinite dilution a fraction a (arbi- 
trarily close to 1) of the fraction /^ stays within a finite 
radius Tq,. It has thus been common practice to call /^ 
the fraction of condensed counterions or Manning frac- 
tion, although limQ,^! r^ = oo. Actually, Rm diverges 
like VR, which follows dir| eGtly from either of the asymp- 
totic boundaries ^ for ^iM3 

Investigating counterion condensation by means of 
computer simulations requires a criterion which identi- 
fies condensed ions. Here we show that the functional 
form of the counterion distribution function suggests a 
simple rule for recognizing the Manning layer. 

If the function P is known, the condensed counterion 
fraction can be characterized in the following "geomet- 
ric" way: Eq @ shows that P viewed as a function of 
ln(r) is merely a shifted tangent- function with its cen- 
ter of symmetry at {ln(_RM); f^}- Since tan"(0) = 0, the 
Manning radius and Manning fraction can be found by 
plotting P as a function of ln(r) and localizing the point 
of inflection. 

This property of P, derived within the framework of 
PB-theory, can in turn be used to define the condensed 
fraction. It provides a suitable way to quantify counter- 
ion condensation beyond the scope of PB-theory, and it 
is exact in the PB-limit without added salt, by construc- 
tion. 

Our counterion condensation criterion can be reformu- 
lated in terms of the counterion density n{r): If P has a 
point of inflection as a function of Inr, dP/d ln(r) must 
have a stationary point there. Using eqs || and |9| it fol- 
lows that r^ n{r) must have a stationary point, which in 
the simple salt-free case is actually a minimum. (In fact, 
for our simulated data we localized the point of inflection 
by fitting a (2,2) Pade approximant in Inr to r"^ n{r) in 
the vicinity of its minimum.) 

It is appropriate to mention bricfiy here three other 
methods which have been used to measure counterion 
condensation and point out their shortcomings. The no- 
tion of a condensed layer closely surrounding the rod 
suggests determining the condensed fraction by simply 
counting the ions within a certain (small) distance of the 
rod, say, a few diameters or one screening length 1/kJIJ 
This amounts to making a prior assumption about the 
Manning radius. Such a procedure is not only arbitrary; 
moreover, the PB Manning radius depends on the poly- 
electrolyte density and diverges like y/R in the dilute 
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FIG. 1. Three predictions for the fraction of condensed 
counterions for a cell with ro = cr and ^ = 2 as a function 
of cell size R. The solid line is the inflection point criterion. 
The dashed line derives the distinction between condensed 
and uncondensed ions from the condition y(r) = 1 while the 
dotted line takes y(r<i] — y{r) — 1. The two arrows mark 
the values for the energy based criteria in the limit 7? — > oo. 
Notice that they do not coincide with the Manning fraction 
1 - 1/e = 1/2. 



limit. If this is not taken into account, the condensed 
fraction is either underestimated (for fixed condensation 
distance) or overestimated (for a distance proportional 
to the screening length 1/k of the counterions, which is 
proportional to R). 

Conversely, one could assume that the condensed frac- 
tion is always given by 1—1/^ and thereby obtain the 
size of 1|he condensed layer, e.g. when salt is added to the 
systemJla Although being exact in the salt-free PB limit, 
this criterion excludes by definition the possibility that 
any effects beyond the mean-field level (like correlations) 
or the presence of salt also modify the fraction of con- 
densed counterions. It also does not predict a crossover 
to a high salt regime where all counterions are condensed 
solely due to the presence of the salt (see section VI). 

Finally, one could be tempted to bring into play the 
electrostatic binding energy and regard all ions within 
a thermal|-diS|tance R^ defined by yiRr) = 1 as 
condensed.llSE3 (Alternatively, one might require a po- 
tential difference of fceT/eo with respect to the rod sur- 



face, i.e., y{ro) — y{R'rj,) = 1.) Within salt-free PB theory, 
however, this is not a suitable criterion since the value of 
the electrostatic potential at the Manning radius, ^{Rm), 
is in no way special (upon dilution it actually diverges 
logarithmically with respect to the boundary aa-|Well as 
the rod surface, as can be derived from cq [k]) E3 As an 
illustration, figure m compares the inflection point rule 
and two energy based criteria with regard to their pre- 
dictions for the condensed fraction: All three methods 
quantify condensation differently, and the energy based 
approaches are density dependent - which on its own is 
not a problem. The unsatisfying aspect is rather that 
the latter do not converge against the Manning frac- 
tion 1 — 1/^ upon dilution. In fact, using eq nu it can 
e.g. be shown that i?T ~ ci i? in the dilute limit, with 
ci ~ 0.25726 being a solution of ci(l — In ci) = e~^/^. 



Observe that i?T scales linearly with R and thus faster 
than the Manning radius; therefore, the thermal radius 
will enclose more than the Manning fraction in the di- 
lute limit. In fact, P{Rt) ~ 1 — 02/^, with the con- 
stant C2 = In ci/(ln ci — 1) R:! 0.57585. Incidentally, this 
implies the relative deviation of this approach from the 
Manning fraction 1 — 1/^ to become small at large £,. 

Let us thus repeat that the inflection point criterion 
employed in the present work has the advantages of (i) 
not fixing by definition the amount of condensed coun- 
terions (/j and Rm can be determined independently 
of each other), (m) reproducing the salt-free PB limit, 
namely P{Rm) = 1 — 1/Ci ^ind (iii) quantifying the break- 
down of the coexistence of condensed and uncondensed 
countcrions in the high salt limit, as will be shown in 
Sec. 0. 

In concluding this section we note that the appearance 
of the logarithm in the inflection point criterion is related 
to ln(r) being the 2D Coulomb potential, i.e., the Green 
function of the cylindrically symmetric Laplacian. In the 
corresponding 3D (spherical) problem of charged colloids 
the Green function 1/r would be Ac appropriate choice 
for plotting the radial coordinate. t3 



IV. COMPARISON OF PB THEORY WITH 
SIMULATIONS: NO ADDED SALT 

In this section we supplement the results of salt-free 
PB theory with computer simulations of a cell-like model, 
with particular emphasis on the role of Manning param- 
eter and valence. Details of the model, the simulations 
and our notation conventions are summarized in the ap- 
pendix. 

Figure shows the counterion distribution functions, 
P{r), for three systems with monovalent countcrions, 
ie/ro = 1, R/ro = 123.8 and ^ G {0.96,1.92,2.88}, i.e., 
counterion condensation is expected to occur for the lat- 
ter two. As suggested in the previous section the func- 
tions are plotted using a logarithmically scaled r-axis. 
Note that in all our PB calculations and simulations the 
distance of closest approach to the rod was tq = ct, where 
a is the small ion diameter, used in the simulations (see 
appendix) . 

For the system with ^ — 0.96 in figure || the agree- 
ment between simulation and PB theory is excellent - 
deviations are almost within the linewidth of the plotted 
curves. For the other two cases the agreement fails quan- 
titatively, but not qualitatively: The shape of the distri- 
bution function remains largely unchanged. The curves 
for P{r) lie above the PB result for all r, indicating a 
stronger condensation than predicted by mean-field the- 
ory. Importantly, for ^ > 1 the simulated curves display 
a point of inflection. As described in the previous section 
this can be used to define a Manning radius and a con- 
densation fraction, which permits one to quantify by how 
much the condensation is stronger. This is summarized 
for a range of densities in table 1, where it can be seen 
that deviations towards higher condensation are stronger 
for dense systems and relax towards the PB prediction 
upon dilution. 

Note that in figure y the measured fraction /^ is larger 







FIG. 2. Simulated counterion distributions (solid lines) 
and PB-results (dotted lines) for a monovalent system with 
ro = la, R/ro ~ 123.8 and (from bottom to top) 
£, G {0.96, 1.92, 2.88}. The f-arrows mark the inflection points 
in the PB-distribution while the ^-arrows mark those points in 
the MD distributions. Note the logarithmically scaled r-axis 
in the present and the following figures. 



than the PB-prediction, while the corresponding Man- 
ning radius is smaller, but since within PB-theory Rm 
increases monotonically with ^, the measured curves can- 
not be modelled by a PB distribution with a somewhat 
larger effective ^. 

As already mentioned, within PB theory the shape of 
the integrated distribution function, P(r), depends on 
the Manning parameter S, and the counterion valence v 
only via the product ^v. We now show that this is a 
property seen on the mean-field level only. In addition 
to the PB result and the simulation of the system with 
^ = 2.88 already shown in flgure 0, the simulated P{r) 
curve for a system with trivalent countcrions and ^ — 
2.88/3 = 0.96 is displayed in flgure 0. The corresponding 
P{r) indicates an even stronger condensation than the 
simulation for ^ = 2.88 and monovalent countcrions. 

There are two principle ways in which PB theory 
can fail: (i) the neglect of excluded volume interactions 
("point-like" ions) and (ii) missing correlations. For the 
first point there is a simple self-consistency check: The 
counterion density is highest at the surface of the inner 
rod and e.g. in the limit i? — > cxd given by eq|lj, but coun- 
tcrions with a finite size might not be able to give rise 
to a density as large as that. This consequently limits 
the range of applicability of the PB approximation to- 
wards not too large ions, not too small cylinders and not 



R/ro 



C = 1.92 

5 = 2.88 
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0.56 
0.78 



0.57 
0.76 



0.56 
0.73 



0.54 
0.70 



0.53 
0.69 



0.51 
0.67 



0.49 
0.67 



TABLE I. Measured condensation fraction /^ for vari- 
ous monovalent systems with ro — (J, which differ in cell 
size R and thus polyelectrolyte density. Within PB theory 
/^ — 1 — 1/^, giving /1.92 ~ 0.479 and /2.8s ~ 0.653 - indepen- 
dently of R. Note that the counterion distribution functions 
for R/ro = 124 are shown in figure H. 




FIG. 3. Comparison of Manning parameter and valence. 
The simulated system from figureH with ^ = 2.88 (sohd curve) 
and its PB-solution (dotted curve) are contrasted with a sim- 
ulation where the Manning parameter is three times smaller 
but the counterions are trivalent (dashed curve). 



too strong electrostatics. In our simulations we are below 
that limit. Excluded vol urt i c interactions - if present - re- 
duce the contact densityBcj and hence could only lower 
the integrated distribution function P(r)] this however 
is not found in the simulations. The observed stronger 
condensation must therefore be attributed to correlations 
neglected in the PB approach. Since an increase in den- 
sity goes along with an increase of correlations, this ex- 
planation seems to be intuitively correct and it is also 
supported by quasi-aoalytical theories, which go beyond 
the mean-field level.ESEl 



V. COUNTERION CONDENSATION IN THE 
PRESENCE OF SALT 

The central question to be discussed in this section is 
the following: Can the concept of counterion condensa- 
tion be extended to the case of added salt? 

Note first of all that the salt corresponds to a new de- 
gree of freedom which comes along with its own length 
scale, namely, a Debye length ^d — {8T:£Bv'^n)~^^^, where 
V is the valence of the (for simplicity) symmetric salt 
and n is its density. It is of central importance how 
this new length relates to the characteristic length Rm 
of the condensation structure: If £d is large compared to 
Rm-, the pure-counterion n{r) structure is preserved; if 
^D is smaller, it dictates the shape of the charge distribu- 
tion function and the condensation structure is no longer 
present. 

Since counterion condensation becomes apparent in 
the behavior of the charge distribution function for R —^ 
oo, one should also investigate this limit in the presence 
of salt. This however is crucially dependent on the chosen 
ensemble^ i.e., whether the limit is performed at constant 
number N of salt molecules or at constant chemical po- 
tential fi. 

In the constant N case €d is proportional to R and 
thus diverges faster than the Manning radius Rm, which 
only scales like \/R. For sufficiently large R the conden- 



sation structure will therefore be visible and the conden- 
sation criterion will be the same as in the salt-free case. 
Conversely, for sufficiently high density or number of salt 
molecules, £d will be smaller than the Manning radius, 
thus modifying the condensation structure. Since the lat- 
ter is a new mechanism for compensating the rod charge, 
it is no longer sensible to use the concept of Manning con- 
densation in this limit. It remains the task of clarifying 
the crossover from counterion condensation to screening, 
which is subject of the following section. 

This line of reasoning needs a little modification if the 
added salt has a higher valence than the counterions, 
since then it will preferably be the salt ions which will 
condense onto the rod. Two cases have to be distin- 
guished: 

1. Already a fraction of the negative salt ions of high- 
est valence could completely neutralize the rod. If 
these ions are taken to be the "true" counterions 
and all the rest (including the "original" counteri- 
ons) is denoted as "salt", one can expect a Man- 
ning limiting behavior typical for the highvalent 
new counterions. 

2. There is not enough salt to completely neutralize 
the rod with the negative salt ions. This is just 
as complicated as the salt-free case with different 
species of counterions and will not be pursued fur- 
ther in this paper. 

Quite differently, in the constant /i case the Debye length 
of the salt will remain finite in the limit R ^ oo and 
consequently smaller than the diverging Manning radius. 
The condensation structure will always be wiped out in 
the infinite dilution limit and it is not possible to produce 
a condensation criterion along the lines of the salt-free 
case. We therefore prefer to work in the constant N 
ensemble. 

From an experimentalists point of view, keeping N or 
/I constant in the limit i? — > cx) corresponds to two com- 
pletely different procedures: In the first case the polyelec- 
trolyte solution is diluted by the addition of pure water. 
In the second case the dilution is done with a salt solu- 
tion of the same ionic strength as the one in which the 
polyelectrolyte originally has been dissolved. One there- 
fore cannot expect these two cases to become equivalent 
in the thermodynamic limit. 

In large systems it is essentially irrelevant whether a 
certain salt concentration is achived by choosing a certain 
number of salt ions or a corresponding chemical potential 
for them. Yet, in all typical systems accessible to com- 
puter simulations the number of ions is still rather small, 
so that the chosen ensemble matters. Since eventually 
a comparison between simulation and theory is headed 
for and since the most straightforward ensemble for sim- 
ulations is the one which conserves particle number, the 
following section derives a PB equation in the presence 
of salt for this case. 



VI. POISSON-BOLTZMANN EQUATION FOR 
CONSTANT NUMBER OF SALT MOLECULES 

Assume that in addition to the monovalent counterions 
of the positively charged rod the cell contains K different 
V : V salts of concentrations fiy with v = I . . . K. The 
overall concentration of negative monovalent ions is thus 
fii + m with m = X/enR'^ where m is the contribution 
due to the counterions of the rod. 

The free energy F = U — TS accounts for the internal 
electrostatic energy U and the translational entropy S of 
the mobile ions in solution. It can be written in terms of 
the electrostatic potential $ and the local ion concentra- 
tions Uv and n_i, of positive and negative ions of valence 
V, respectively. Within mean-field theory, F is given by 



F ■ 



d\ 



i{v*)' 



K 



k^T y^ Tiy In -^ 

■^— ' 71... 



„^0 



(14) 



where fi_^ — n^ (for v = 2 . . . K) and fi_i = ni + to de- 
note the concentrations of the negatively charged mobile 
ions. 

As discussed in the above section, we are interested in 
the constant N ensemble, i.e., the case that for each ionic 
species the number of ions within the cell of volume V is 
conserved. The local equilibrium concentrations n„ have 
thus to be derived under the constraints 



{riv) = — d^r n,„ 
V 



(15) 



The usual variation of F results then in the Boltzmann 
distributions for the local concentrations 



TlvC 



-vy-n^ 



(16) 



where the chemical potentials /i„ = ln(e^^^) ensure par- 
ticle conservation. 

Again we consider the rod sufficiently long that we can 
neglect end effects. Then, the electrostatic potential $ 
and the local ion concentrations riy depend only on the 
radial distance r to the rod axis. 

Insertion of the local concentrations riy into the cylin- 
drically symmetric Poisson equation e ($" -I- 4>'/r) = 
— ^^ vn^ leads to the Poisson-Boltzmann equation 



-vy 



I K 
— — —AttIb y VHy- ^ 



v=-K 



(17) 



This equation has to be solved subject to the boundary 
conditions (0), i.e. the same as for the salt-free case. 



Numerical solutions of eq 17 can be found employing a 
Newton-Raphson iteration scheme in which the chemical 
potentials /x„ = ln(e~''^) are updated after each iteration 
step. Once a solution y{r) is found, the integrated charge 
distribution function of the mobile ions 
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P{r) = 1 / dr' 27rr' 
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(18) 



can simply be calculated by P{r) = 1 + ry'{r)/2^ which 
follows from inserting the Poisson equation into eq Vn 
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FIG. 4. PB results (dotted curves) for the integrated 
charge distribution function P{r) for a system characterized 
by ro = cr, R/ro = 61.9, ^ = 2.1 and A = 0.96 e/ro. Note that 
the number of rod counterions corresponding to rod length 
L = 250ro is M = 240, while the number of 1:1 salt coun- 
terions (and co-ions) per length L is, from bottom to top, 
TV = 0, 104, 800, 3070 and 15 000. The bold solid curve shows 
the locus of inflection points, i.e., the union of all inflection 
points of the functions P(r). The f-arrow marks the loca- 
tion of the salt-free Manning inflection point and the |-arrow 
shows where it joins one of the new salt inflection points. 
The branch of the locus between these two arrows indicates 
the range in which the concept of Manning condensation is 
meaningful. Observe that the functions are convex-up only 
within the grey-shaded region. 



and carrying out the integration with consideration of the 
boundary conditions (^). Since P{r) is a measure of the 
fraction of the overall electrolyte charge found within a 
cylinder of radius r, we must have P{ro) = and P{R) — 
1. Note that eq |l^ is a natural generalization of the 
distribution function from eq 0, but its interpretation 
as an integrated probability distribution (or fraction of 
counterions) is only valid in the salt-free case. 

To investigate the condensation criterion in the pres- 
ence of monovalent salt we have calculated the mean- 
field potential y{r) solving the PB equation for a system 
characterized by tq — a, R/ro — 61.9, A/e — 0.96/ro, 
^ = 2.1 and a variable number of salt molecules. To fa- 
cilitate the comparison with computer simulations in the 
next section we write TV for the number of monovalent 
salt molecules associated with a rod segment of length 
L = 250 rg. The corresponding cell volume that contains 
the mobile ions is then V = LttR'^ and the Debye length 
is .^D = {8TriBN/Vy^^^. Note that the hne charge den- 
sity A = 0.96 e/ro implies a number M = 240 positive 
charges found on the rod segment of length L = 250 tq, 
and that the investigated N ranges from to 15 000. 

From the numerical solutions of y{r) we have deter- 
mined all inflection points of P{r) plotted against Inr. 
These inflection points are solutions of the equation 
d^P(r)/d(lnr)^ = 0. In figure^ we present the inflec- 
tion points (bold solid curve) starting from JV = up to 
N = 3070. For larger values of N no further inflection 
points arc found. We also show the integrated charge dis- 
tributions (dotted) for A^ = 0, 104, 800, 3070 and 15 000, 
corresponding to Debye lengths oi £u/ro — oo, 22.9, 8.3, 







FIG. 5. Simulated curves P{r) for the systems with A*' = 0, 
104 and 3070 from figure (solid lines) and the corresponding 
PB results (dotted curves). 



4.2 and 1.9, respectively. 

The location of the inflection point for A^ = coincides 
with i?M; thus indicating a fraction of condensed counte- 
rions of P{Rm) = 1 — 1/^- Increasing N by adding salt 
shifts the inflection point to smaller values of r. That 
is, the layer of condensed counterions contracts, which 
is in accord with other condensation criteria mentioned 
Importantly, the amount of condensed 



in section [II 



counterions is only marginally increased in the presence 
of monovalent salt. From a certain N on (in figure H 
we find A^ = 104) two more inflection points appear in 
the high r/rg region. This happens typically for a corre- 
sponding Debye length being of the order of the cell size 
itself, indicating the appearance of a characteristic, salt 
induced, change in the convexity of P (as a function of 
Inr). 

Upon a further increase in N one of the two new inflec- 
tion points shifts towards smaller r/ro values, finally fus- 
ing with the Manning inflection point. Roughly speaking, 
we find the inflection points to vanish if the Debye length 
characterizing the salt content becomes smaller than the 
radius of the condensed layer. This suggests a break- 
down of the need to distinguish between condensed and 
uncondensed counterions once the typical salt screening 
length interferes with the size of the condensed counte- 
rion layer. Indeed, for a very high salt content, where the 
Debye length is much smaller than the radius of the rod, 
the solution of the PB equation would be the one of a 
charged plane and one may consider all excess counteri- 
ons being condensed no matter what the charge density 
of the rod is. 



VII. COMPARISON OF PB THEORY WITH 
SIMULATIONS: ADDED SALT 

In this section we again compare the numerical results 
of the PB equation - this time in the presence of salt 
- with computer simulations. We reinvestigate the sys- 
tems in figure with ^ — 2.1, monovalent counterions 
and number of salt molecules A^ = 0, 104, and 3070 with 




FIG. 6. The integrated charge distribution function P{r) 
for a system characterized by ro — a, R/ro — 15.5, ^ = 4, 
A/e = 0.96/ro, A^ = 1000 molecules of a divalent salt and 
M = 60 monovalent counterions corresponding to a rod seg- 
ment of length L = 62.5ro. The simulation (solid curve) shows 
a pronounced overcharging-effect, in contrast to PB-theory 
(dotted curve). 
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FIG. 7. The densities n_2(r) (solid line) and n+2{r) (dot- 
ted line) of negative and positive salt ions, respectively, for 
the same system presented in figure M. The inlay shows the 
electrostatic potential y(r). 



respect to a rod segment of length L = 250ro. The re- 
sults of the computer simulations and the corresponding 
mean-field calculations are presented in figure @. Like in 
the salt-free case the computer simulations show a more 
pronounced condensation effect towards the rod, which 
we again attribute to ion-ion correlation effects. Still, 
the shape of the distribution functions remains qualita- 
tively the same. Note in particular that the appearance 
and disappearance of two points of infiection at iV = 104 
and N — 3070 respectively, which leads to extremely 
small curvatures in the PB distribution functions, also 
leads to very straight regions in the measured distribution 
functions. The crossover from Manning condensation to 
screening, as described within PB theory, can thus be 
expected to be essentially correct. 

The PB approach fails to describe the physical situ- 
ation if one or more of the following points apply: (i) 



the electrostatic interactions are strong, (ii) the coun- 
terions are multivalent and (Hi) the density is high. A 
simulation under such conditions can be inspected in fig- 
ure |g. In this system we have tq = a, R/tq « 15.5, 
A ~ 0.96e/ro, ^ = 4, 60 monovalent counterions, and 
1000 molecules of a 2:2 salt (giving a Debye length of 
roughly 0.33 tq, i.e., smaller than the ion diameter). Here 
P{r) overshoots unity, showing a charge reversal of the 
rod at distances around r « 1.5 rg, while the simple PB 
prediction is clearly qualitatively off. This phenomenon 
is usually referred to as overcharging and has been pre- 
dicted for the primitioie cell model first from hypernet- 
ted chain calculation;^ and later by a modified Poisson- 
Boltzmann approach.Ej 

Since P{R) = 1, the overshooting above 1 at small 
distances implies the existence of a range of r-values at 
which the mobile ion system is locally positively charged 
(i.e., with the same charge as the rod), such that P{r 
can eventually decay to 1. This is seen in figure 
which shows that n+2{r) > n^2{f') at r « 2ro. Since 
P(1.5) « 1.45, the rod and its innermost layer of con- 
densed ions could be viewed as an effective rod of radius 
1.5 ro which is negatively charged with Manning param- 
eter ^ = 1.8. Since this value is again larger than 1, it 
entails ion condensation (this time of the positive ions). 
In fact, it even leads to a second overcharging, as can 
clearly be seen in figure o, where P{r) - in decaying from 
1.45 - overshoots the value of 1 again. Overcharging can 
thus give rise to layering; in the presented example no 
less than three layers can clearly be made out. These 
local charge oscillations also reflect themselves in oscil- 
lations of the electrostatic potential, as demonstrated in 
the inset of figure M. Note that these oscillating poten- 
tials will also have pronounced effects on the interaction 
between such rigid polyelectrolytes. 



Upon addition of monovalent salt we found PB the- 
ory to predict counterion condensation within a some- 
what smaller region around the charged rod. This is 
in accord with|-0|ljlifpstudicsE3EZl using different conden- 
sation critcria.ll3l21]'E3 Importantly, the fraction of con- 
densed counterions did not exhibit a strong salt depen- 
dence thus supporting a principal idea in Manning con- 
densation: The number of condensed counterions does 
not depend on salt. Yet, this observation must break- 
down in the high salt limit, where the screening length of 
the salt becomes of the order of the size of the condensed 
counterion layer. Indeed, our new condensation criterion 
naturally predicts this to happen. 

All these behaviors were well reproduced by our com- 
puter simulations. In fact, the agreement between simu- 
lation and the corresponding mean-field level calculation 
is remarkably good for systems with a Manning param- 
eter ^ < 1. Upon increasing the counterion density, va- 
lence, or the Manning parameter, the simulations predict 
consistently a somewhat stronger condensation. We have 
argued that this finding is due to ionic correlations not 
present in PB theory. To test this assumption, we have 
also performed a simulation of a system in a highly con- 
centrated divalent salt environment. Here we clearly saw 
the phenomenon of overcharging, which a corresponding 
PB calculation was unable to reproduce. 

We note finally that our MD simulations suggest the 
usability of the counterion condensation criterion beyond 
the PB cell model approximation. In fact, we have ob- 
served inflection points in the integrated probability dis- 
tribution functions in systems very distinct from the rigid 
rods employed in the present study, namely for flexible 
polyelectrolytes in the presence of multivalent counteri- 
ons and for flexible polyelectrolytes in poor solvents with 
monovalent counterions.tj 



CONCLUSIONS 
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We have revisited counterion condensation with and 
without added salt for a solution of rigid polyelectrolytes 
within the cell model approximation. It was confirmed 
that on the level of PB theory with no added salt a simple 
geometric method locates the condensation radius Rm as 
well as the fraction of condensed counterions. This ge- 
ometric method consists in finding the infiection point 
of the integrated probability distribution P plotted as 
a function of Inr. Without added salt, the locations of 
the inflection point and the Manning radius Rm, where 
a fraction of 1 — l/<^ counterions are condensed, are iden- 
tical. 

A key point in the present work was to extend the 
inflection point based counterion condensation criterion 
to the case of added salt and to compare its implications 
as predicted by PB theory and by computer simulations. 
Our motivation for introducing this new condensation 
criterion was (i) to avoid fixing by definition the amount 
of condensed counterions, (ii) to reproduce the salt-free 
PB limit, namely P{Rm) = 1 — l/C? and (Hi) to predict a 
counterion condensation breakdown in the high salt limit. 
We are not aware of any other condensation criterion that 
fulfills all these requirements at the same time. 
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APPENDIX: DETAILS OF THE SIMULATION 

The system used to study counterion condensation 
consists of a cubic simulation box of length Lb, a charged 
rod along the main diagonal, the appropriate amount of 
counterions necessary for electric neutrality and possi- 
bly some additional salt. Upon switching on 3D-periodic 
boundary conditions a triangular array of infinitely long 
charged rods is modelled. Note that in PB theory we 
approximate the corresponding Wigner-Seitz-cell by a 
cylindrically symmetric unit cell of the same volume (im- 
plying R = Lb/v^W^) thus rendering the PB equation 
one-dimensional. 

Apart from electrostatic interactions all ions are sub- 
ject to a purely repulsive Lennard- Jones potential, giving 



an excluded volume and a corresponding ion-diameter a. 
Further, the rod is modelled as an immobile string of such 
spheres, having the separation 1.042cr. The distance of 
closest approach to the cylinder, i.e. tq, turns out to be 
essentially 1 a. 

The electrostatic interactions in this periodic boundapj; 
geometry were computed wijtli the help of P'^M routineg23 
and a Langevin thermostatic combined with a velocity- 
Verlet-integratoiEj (with timestep 0.0125 in LJ units) 
was implemented to drive the system into the canonical 
state. The number of MD steps varied between 8 x 10^ 
(for the systems from figure |) up to 6.4 x 10^ (for the 
system from figures ^ and |7) and the saturation of the 
electrostatic energy was used to test for equilibration. A 
more detailed description of our simulatiop-piethod will 
be presented in a forthcoming publication.c3 
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